clc;
clear all;
close all;

%% 测试参数设置
sizes = [100,200,300,400,500,600]; % 测试不同矩阵规模
r = 50; % 固定秩
num_sizes = length(sizes);

% 预分配结果存储
times_pinv = zeros(num_sizes, 1);
times_qril = zeros(num_sizes, 1);
times_direct = zeros(num_sizes, 1);
times_schmidt = zeros(num_sizes, 1);

errors_pinv = zeros(num_sizes, 1);
errors_qril = zeros(num_sizes, 1);
errors_direct = zeros(num_sizes, 1);
errors_schmidt = zeros(num_sizes, 1);

%% 主测试循环
for i = 1:num_sizes
    n = sizes(i);
    m = n; % 方阵测试
    fprintf('\n===== 测试矩阵大小: %d x %d =====\n', m, n);
    
    % 生成固定秩的随机矩阵
    A = randn(m, r) * randn(r, n); 
    T = eye(m);
    
    % 清除可能的冲突变量
    clear QRIL_inverse
    
    %% 算法执行与计时
    % MATLAB内置pinv
    tic;
    A_pinv = pinv(A);
    times_pinv(i) = toc;
    
    % 自定义QRIL算法
    tic;
    QRIL_inv = QRIL_inverse.QR_schmidt(A, T);
    times_qril(i) = toc;
    
    % directly_MLE_iter
    tic;
    Directly_MLE_inverse = QRIL_inverse.directly_MLE_iter(A, T);
    times_direct(i) = toc;
    
    % Schmidt_MLE_solver
    tic;
    Schmidt_MLE_inverse = QRIL_inverse.Schmidt_MLE_solver(A,T);
    times_schmidt(i) = toc;
    
    %% 计算误差
    errors_pinv(i) = norm(A*A_pinv - T, 'fro');
    errors_qril(i) = norm(A*QRIL_inv - T, 'fro');
    errors_direct(i) = norm(A*Directly_MLE_inverse - T, 'fro');
    errors_schmidt(i) = norm(A*Schmidt_MLE_inverse - T, 'fro');
    
    %% 显示当前结果
    fprintf('计算时间: pinv=%.4fs, QRIL=%.4fs, MLE_iter=%.4fs, Schmidt_MLE=%.4fs\n',...
        times_pinv(i), times_qril(i), times_direct(i), times_schmidt(i));
    fprintf('误差范数: pinv=%.4e, QRIL=%.4e, MLE_iter=%.4e, Schmidt_MLE=%.4e\n',...
        errors_pinv(i), errors_qril(i), errors_direct(i), errors_schmidt(i));
end

%% 绘制性能曲线
figure('Position', [100, 100, 1200, 800]);

% 时间性能曲线
subplot(2, 1, 1);
semilogy(sizes, times_pinv, 'b-o', 'LineWidth', 2, 'MarkerSize', 8);
hold on;
semilogy(sizes, times_qril, 'r-s', 'LineWidth', 2, 'MarkerSize', 8);
semilogy(sizes, times_direct, 'g-d', 'LineWidth', 2, 'MarkerSize', 8);
semilogy(sizes, times_schmidt, 'm-^', 'LineWidth', 2, 'MarkerSize', 8);
grid on;
title('求逆算法时间性能比较', 'FontSize', 14);
xlabel('矩阵规模', 'FontSize', 12);
ylabel('计算时间 (s, log scale)', 'FontSize', 12);
legend('MATLAB pinv', 'QRIL', 'MLE\_iter', 'Schmidt\_MLE', 'Location', 'northwest');
set(gca, 'FontSize', 11);

% 误差性能曲线
subplot(2, 1, 2);
loglog(sizes, errors_pinv, 'b-o', 'LineWidth', 2, 'MarkerSize', 8);
hold on;
loglog(sizes, errors_qril, 'r-s', 'LineWidth', 2, 'MarkerSize', 8);
loglog(sizes, errors_direct, 'g-d', 'LineWidth', 2, 'MarkerSize', 8);
loglog(sizes, errors_schmidt, 'm-^', 'LineWidth', 2, 'MarkerSize', 8);
grid on;
title('求逆算法精度比较', 'FontSize', 14);
xlabel('矩阵规模', 'FontSize', 12);
ylabel('Frobenius 误差范数 (log scale)', 'FontSize', 12);

legend('MATLAB pinv', 'QRIL', 'MLE\_iter', 'Schmidt\_MLE', 'Location', 'northwest');
set(gca, 'FontSize', 11);
set(gca, 'YScale', 'log')  % 关键设置：y轴对数坐标系
%% 保存结果
save('inverse_perf_results.mat', 'sizes', 'times_pinv', 'times_qril', 'times_direct', 'times_schmidt',...
    'errors_pinv', 'errors_qril', 'errors_direct', 'errors_schmidt');

%% 显示结果表格
fprintf('\n===== 综合性能结果 =====\n');
fprintf('矩阵规模\t pinv时间\t QRIL时间\t MLE时间\t Schmidt时间\t pinv误差\t QRIL误差\t MLE误差\t Schmidt误差\n');
for i = 1:num_sizes
    fprintf('%dx%d\t\t%.4f\t\t%.4f\t\t%.4f\t\t%.4f\t\t%.2e\t%.2e\t%.2e\t%.2e\n',...
        sizes(i), sizes(i),...
        times_pinv(i), times_qril(i), times_direct(i), times_schmidt(i),...
        errors_pinv(i), errors_qril(i), errors_direct(i), errors_schmidt(i));
end